Simulating adiabatic evolution of gapped spin systems 
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We show that adiabatic evolution of a low-dimensional lattice of quantum spins with a spectral gap 
can be simulated efficiently. In particular, we show that as long as the spectral gap /S.E between the 
ground state and the first excited state is any constant independent of n, the total number of spins, 
then the ground-state expectation values of local operators, such as correlation functions, can be 
computed using polynomial space and time resources. Our results also imply that the local ground- 
state properties of any two spin models in the same quantum phase can be efficiently obtained from 
each other. A consequence of these results is that adiabatic quantum algorithms can be simulated 
efficiently if the spectral gap doesn't scale with n. The simulation method we describe takes place in 
the Heisenberg picture and does not make use of the finitely correlated state/matrix product state 
formalism. 

PACS numbers: 03.67.Lx, 75.10.Pq, 75.40.Mg 



I. INTRODUCTION 

The low-temperature physics of lattices of interact- 
ing quantum spins is typically very complex. The com- 
putational cost of even approximating basic properties, 
such as the ground-state energy eigenvalue, of these sys- 
tems is prohibitive. Indeed, for 2D lattices of interacting 
spins, the task of computing an approximation to the 
ground-state energy eigenvalue correct to within some 
polynomial confidence interval is fantastically difficult — 
this problem is complete for the complexity class QMA, 
which is the quantum version of the complexity class N P 
flii. 

It might therefore seem that the computational task of 
approximating the low-temperature behaviour of inter- 
acting quantum spins is entirely hopeless. However, for 
physically realistic models, this is not the case in practice. 
Many algorithms have been developed which appear to 
provide efficient approximations to a wide variety of local 
properties of physically realistic systems, such as correla- 
tors, at low temperature. Perhaps the most successful of 
these methods has been the family of algorithms based on 
the density matrix renormalisation group (DMRG) (See 
[3| and references therein for a review of the DMRG and 
description of extensions.) 

The DMRG is a remarkably flexible and adaptable 
algorithm, admitting a slew of generalisations. Appli- 
cations include: simulating dynamics 0, Q, dissipative 
systems [tI, Isj , disordered systems [oj , and higher dimen- 
sional lattices J^]. At least part of the flexibility of the 
DMRG is due to the fact that it is equivalent to a vari- 
ational minimisation over the space of finitely correlated 
states (FCS) Hence, the methodology of the DMRG 
can be adapted to any situation where the principle ob- 
ject of study, be it an eigenstate or propagator, can be 
approximated using a FCS vector on Hilbert space. An 
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alternative to methods based on variations over FCS has 
been recently proposed which appears to offer spectacu- 
lar computational speedups over the DMRG and relatives 

In practice it appears that the DMRG and related al- 
gorithms can efficiently obtain arbitrarily accurate ap- 
proximations to the local ground-state properties of a 
ID collection of interacting quantum spins. However, at 
the current time, there is no satisfactory understanding 
of the correctness (i.e. will the DMRG always return a 
faithful approximation to the ground state and not some 
other eigenstate) and the complexity (i.e., assuming cor- 
rectness, how much computational resources are required 
to obtain a good approximation to a ground state) of the 
DMRG. 

The correctness of the DMRG is far from obvious. 
This is because the ground-state approximation obtained 
by the DMRG cannot be certified; the DMRG only re- 
turns an approximate ground-state eigenvector and can- 
not guarantee that this vector is close to the true ground 
state. It is therefore extremely desirable to determine 
a priori the class of systems for which the DMRG and 
relatives provably return faithful approximations to the 
ground state. The complexity [40] of the DMRG is also 
difficult to ascertain. Assuming we could prove correct- 
ness of the DMRG for a class of realistic physical systems, 
the actual complexity of the DMRG depends subtly on 
many detailed properties of the system, such as geomet- 
ric entropy of the ground state, and nonconvexity of the 
objective function which is minimised. 

Recently this situation is changing [l^, [3, ■ In [T^ 
an analysis of the resource scaling of a DMRG-like al- 
gorithm to obtain approximations to the ground states 
of ID gapped local models was undertaken. This pa- 
per provides the first general subexponential estimate for 
the time and space resource requirements of any prov- 
ably correct method to compute approximations to the 
ground states of gapped models; it was found that if the 
model is gapped then resources scaling as rt'^'°s", with c 
some constant, are sufficient to obtain and store a com- 
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putational representation of the ground state of a gapped 
local model E|. In [l^ it was shown that the ground 
state of some |42| critical ID spin models can be stored 
efficiently. Unfortunately, there is currently no theoreti- 
cal argument which implies that these approximations to 
the ground states can be obtained efficiently. Indeed, the 
results of this paper imply that if such approximations 
are obtained via adiabatic continuation then exponential 
computational resources may be required to obtain them. 
(Note, however, that we can say nothing about the other 
methods to obtain such FCS approximations.) Finally, 
in it was shown that an approximation to the prop- 
agator for a ID lattice of quantum spins can be obtained 
and stored (as a FCS vector) using polynomial resources 
in n and the error e and exponential resources in the time 
\t\. (It is straightforward to extend the argument of 
to show an analogous result in 2D.) 

There is at least one solid reason why we believe that 
DMRG-like methods ought to provide a computationally 
efficient recipe to compute approximations to the ground 
states of gapped systems. Namely, we know that the 
ground-state correlation functions for any gapped sys- 
tem are clustering or rapidly decaying with separation 
[13) dj 13 ■ This result, which is the natural analogue 
of Fredenhagen's proof [131 of clustering for relativis- 
tic quantum field theories, is especially impressive given 
that it applies to an extremely wide class of quantum 
lattice systems in low dimensions. As a consequence of 
clustering results we conclude that gapped spin systems 
are essentially free — an intuition which is persuasively 
backed up by classical renormalisation-group style argu- 
mentation — and thus can be modelled as noninteracting 
effective spins, which can be simulated easily. 

Another way of arriving at this conclusion is to think of 
correlations as roughly "measuring" the degree of quan- 
tum correlations in the ground state. Since the amount 
of quantum correlations in a quantum state limits the ex- 
tent to which a state can be approximated by a FCS [13] , 
we are strongly encouraged to think that the clustering 
results may actually imply that DMRG-like algorithms 
may converge rapidly for at least some realistic gapped 
systems. 

Unfortunately, knowing that the correlations decay is 
not enough information to infer that the eigenstates are 
finitely correlated. To understand this simply consider a 
generic quantum state [21] which is a quantum state cho- 
sen uniformly from Haar measure induced on state-space. 
A generic quantum state exhibits rapidly decaying corre- 
lations (indeed, all m-point correlation functions are es- 
sentially zero for m < ^) yet such a state is extremely en- 
tangled and cannot be efficiently represented as a finitely 
correlated state. Nevertheless, it might be argued that 
the results of [13, El, [13 avoid this counterexample be- 
cause they prove something stronger, namely exponential 
clustering, which says that the reduced density operator 
PAB of the ground state for two arbitrarily large sepa- 
rated regions A and B is indistinguishable from a product 
Pa ® Pb when it is used to compute expectations for prod- 



uct observables MaMb- Interestingly, a naive attempt 
to exploit this exponential clustering runs into problems. 
The reason is that there exist highly entangled states 
aAB, called data-hiding states, which exhibit precisely 
these properties {2^ ]. Thus, to prove that the ground 
state of a gapped local hamiltonian is well-approximated 
by finitely correlated state with polynomial resources we 
appear to need more information than that given by cor- 
relation functions. 

Despite some recent progress a solution to the fun- 
damental problem, namely, to prove correctness of any 
algorithm which obtains approximations to local ground- 
state properties for gapped ID models and to further pro- 
vide a polynomial theoretical worst-case estimate on the 
resource requirements such an algorithm, still seems far 
away. Let us summarise the various approaches to find- 
ing approximations to the ground state of a spin model 
and the theoretical obstructions encountered in each of 
these approaches. 

There are at least four ways to obtain an approxi- 
mation to the ground state of a quantum system: (i) 
variation over a class of ansatz ground states; (ii) simu- 
lation of the thermalisation process via imaginary time 
evolution or similar; (iii) approximation of the convex 
set of reduced density operators of translation-invariant 
quantum states; and (iv) adiabatic continuation from the 
ground states of classical spin models. The DMRG is an 
example of the first method, namely it is a variation over 
the class of FCS with fixed auxiliary dimension. Unfor- 
tunately this variation is, in general, nonconvex and it 
has been recently discovered [2^ that hard instances for 
a closely related variation problem can be constructed. 
Thus it seems likely that the DMRG is not correct in 
general. The second approach, namely imaginary time 
evolution, suffers from the shortcoming that an initial 
guess for the ground state ]ri) which has nontrivial 
overlap with the actual ground state is required. If such 
an initial guess is unavailable then the storage require- 
ments of the imaginary time evolution approach could 
be, in the worst case, exponential It seems plausible 
that obtaining such a guess could be as hard as solv- 
ing the original problem. The third method requires an 
exponentially good characterisation of the convex set of 
reduced density operators of translation-invariant quan- 
tum states in order to obtain 0(1) estimates for local 
operators. The final method, which is the focus of this 
paper, suffers from the limitation that it is not known 
if the ground state of an arbitrary gapped spin model 
can be obtained via adiabatic continuation from a clas- 
sical model without encountering a quantum phase tran- 
sition. However, it has been recently proved [2j, [2I, j^^] 
that in the neighbourhood of a classical spin model adia- 
batic continuation will work. Thus, using this approach, 
we are able to provide the first polynomial estimates on 
the resource requirements of a correct method to obtain 
a representation of the ground state of at least a subclass 
of gapped models. 

There is an intimate connection between simulating 
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adiabatic continuation for quantum lattice models and 
simulating quantum computations [3, [l^l • Namely, if adi- 
abatic evolution for an arbitrary 2D lattice model with 
a gap that scales as an inverse polynomial of the system 
size could be simulated efficiently on a classical computer 
then [3] BQPCP, thus obviating the need to design and 
engineer a quantum computer in the first place! Natu- 
rally, our results are nowhere near strong enough to show 
the complexity class inclusion BQPCP, but they do have 
implications for error correction methods for adiabatic 
quantum algorithms. 

A complete theory of quantum error correction for adi- 
abatic quantum algorithms [28| is still being developed. 
For example, for general thermalisation decoherence, we 
really have no idea how to calculate a fault-tolerance 
threshold for adiabatic quantum algorithms (see [1^ and 
[30t for a discussion of quantum error correction and 
fault tolerance.). Presumably a general quantum error- 
correcting code for a quantum adiabatic algorithm would 
involve encoding the adiabatic evolution in a larger sys- 
tem such that the minimum spectral gap encountered 
along the evolution was larger This would mean 

it would cost the environment more energy/unit time to 
induce a transition from the ground state during the evo- 
lution (an "error"). It is natural to assume that the gap 
could be boosted to a large constant, independent of the 
number n of spins, with a polynomial increase in size. 
Our results show that if this were possible then we could 
simulate adiabatic quantum algorithms efficiently on a 
classical computer! Thus, conditioned on the strict com- 
plexity class containment P C BQP, we obtain a bound 
on how large the gap could be boosted by encoding for 
adiabatic quantum algorithms. 

The method we develop in this paper is very closely 
related to the method studied in [31]. In 31 1 the au- 
thors investigate the evolution of local operators under a 
quasi-adiabatic change in a local hamiltonian. As long as 
the hamiltonian has a spectral gap throughout the evo- 
lution, it was found that local operators remained local 
and thus it was possible to say that local gauge invariance 
remains when two hamiltonians are in the same phase. 
Our task is similar: we wish to understand the expec- 
tation values of local operators in the ground state of a 
system that has undergone adiabatic evolution. We wish 
to show that the computation of such expectation values 
can be done efficiently on a classical computer as long as 
the smallest gap encountered during the adiabatic evolu- 
tion is 0(1). While this calculation can be treated using 
quasi-adiabatic evolution and the methods developed in 
[3l| to study such evolutions, we prefer to study exact 
adiabatic evolution. We do this primarily in anticipation 
of the application of these results to studying entropy- 
area laws for systems in the same phase. 

We provide an efficient computational method to com- 
pute the expectation values of local operators in the 
ground states of hamiltonians undergoing exact adiabatic 
evolution, a method which works equally well for hamil- 
tonians with spatially varying interactions. Our method 



does not make use of the PCS formalism. Rather, we 
develop our simulation method in the Heisenberg pic- 
ture, where locality is manifest. Indeed, if we were to 
make use of state representations in the Schrodinger pic- 
ture, i.e. the 2D PCS formahsm (PEPS), we would be 
unable to apply our results because even if we could con- 
struct PEPS approximations to the adiabatically contin- 
ued ground state it is currently unknown how to effi- 
ciently extract expectation values of local operators from 
the PEPS representation. We sidestep this issue by pro- 
viding a ground-state certificate in the form of a spec- 
ification of a local hamiltonian which can be efficiently 
numerically simulated in the Heisenberg picture to ex- 
tract local expectation values. 

The outline of this paper is as follows. We begin in 
^Ullbv introducing the class of local hamiltonians we con- 
sider and stating the problem we wish to solve. In ijllll we 
then show how adiabatic evolution for quantum lattices 
of spins can be described by unitary dynamics of an ef- 
fective local hamiltonian. We use this effective dynamics 
in ijlVl to construct an approximate local dynamics which 
can then be used to efficiently extract local properties of 
the adiabatically continued ground state. We conclude 
with some discussion of our results in ^V] We detail some 
simple properties of compactly supported C°° functions 
in Appendix [XI 



II. FORMULATION 

In this section we introduce the Hilbert space and op- 
erator algebras for the systems we consider. We define 
what we mean by strictly local and approximately local 
hamiltonians. Finally, we specify the computational task 
that will occupy us for the rest of this paper. 

We consider quantum systems defined on a set of ver- 
tices V with a finite dimensional Hilbert space Hx at- 
tached to each vertex x £ V. We always assume that 
V is finite. (There are some minor theoretical obstruc- 
tions which currently preclude a simple extension of our 
results to infinite lattices; we'll discuss this in a further 
paper.) For X C V, the Hilbert space associated to X 
is the tensor product Hx — ^xex '^^i and the algebra 
of observables on X is denoted by Ax = B{Hx), where 
B{7ix) denotes the C*-algebra of bounded operators on 
"Hx with norm 



1^11= sup ||A|V.>||, 
WesCHx) 



(1) 



and S{TLx) is the state space for Tix- We assume that V 
is equipped with a metric d. In the most common cases 
V is the vertex set of a graph, and the metric is given 
by the graph distance, d{x,y), which is the length of the 
shortest path of edges connecting x and y in the graph. 
Finally, by tensoring with the unit operators on Y \ X , 
we consider Ax as a subalgebra of Ay, whenever X C Y. 

We will, for the sake of clarity, introduce and describe 
our results for a collection of n distinguishable spin-i 
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particles. Thus, the Hilbert space Ti for our system is 
given by ?i = . We now fix the metric for our 

vertex set V to be that of a low-dimensional periodic 
lattice L oi n — vertices, where m G N and rj is the 
dimension. Because the case 77 = 2 is the only really 
nontrivial case that interests us, we fix 77 = 2 fi'om now 
on. We refer to vertices as sites and identify each site v 
with its coordinates j = {jx,jy)- Because the lattice is 
periodic we identify coordinates: {jx — m) = {j^ = 0) 
and {jy = m) = {jy = 0). It is entirely straightforward 
to generalise our results to higher-dimensional lattices, 
higher dimensional spins, and to more general lattices. 

We consider a distinguished basis, the standard prod- 
uct basis, for Hv given by |z) = {8)™^o 0j"=o l^0xj„))> 
Zj G Z/2Z. We'll also have occasion to refer to a cer- 
tain orthonormal basis for Ay- we denote by cr° = 
(8)7'=o ^J'=o ' ' "^j ^ Z/4Z, the standard opera- 

tor basis, where ( 1 ? ), = ( ? J ), = ( V ) > 

and = ( J _?i ) , are the Pauli sigma matrices. 

We define the support supp(M ) C F of an operator 
M e Av to be the smallest subset A C 1^ such that 
M e ^Ia, i-e., the smallest subset upon which M acts 
nontrivially. Let M C L and N C L. We define the 
sumset M + N C L of M and N by M + N = {x + 
y I X G M, y G TV} where the addition operation x -f y is 
inherited from the standard addition on L = (Z/m'Z) x 
{^jmX). This operation is the natural generalisation of 
the convolution operation on the real numbers to the 
finite group L. (It is fairly straightforward to generalise 
these operations to more general graphs.) 

We now introduce the family i?(s) of parameter- 
dependent hamiltonians we are going to focus on. To 
define our family we'll initially fix some parameter- 
dependent interaction term /i(s) G Av which has 
bounded norm 0: ||/i(s)|| < 0(1). We think of 
/i(s) as being "centred" on site 0, i.e. we demand that 
G supp(/i(s)). Our family His) of quantum systems is 
then defined by 

^(^) = - E'^j(^)' (2) 

where T^x (respectively, '?^) is the unit translation oper- 
ator which translates the subsystems one site across in 
the X (respectively, y) direction, eg., 

(m—lm—l \ m— Im— 1 

and /ij(s) = T^^Tj-ihis))). While the hamiho- 
nian H{s) generated by this construction is translation- 
invariant, none of our subsequent calculations depend on 
this fact in any serious way. Hence the results of this pa- 
per apply equally to hamiltonians with spatially varying 
interactions. 

We are going to make three simplifying assumptions 
about our hamiltonian H{s). The first is that H{s) is as- 
sumed to be strictly local which means that | supp(/i(s)) 



is an 0(1) constant. The second assumption we make is 
that the interaction h{s) that generates H{s) can be writ- 
ten as h(s) = ft-o + sh' , where ho and h' are two operators 
with 0(1) norm. The final assumption is that the ground 
state is unique and the spectral gap AE{s) between the 
ground- and first-excited states for H{s) satisfies the in- 
equality AE{s) > A, Vs G [0,1], where AE{s) is an 
0(1) constant. Note that the first two assumptions can 
be lifted with a little extra work, however, the assump- 
tion that the gap AE{s) is an 0(1) constant cannot be 
relaxed: the simulation algorithm we present scales ex- 
ponentially with AE(s). 

We will also have occasion to discuss approximately lo- 
cal hamiltonians. Such hamiltonians are obtained in the 
same way as in ([2]), that is, we fix some initial interac- 
tion term k{s) which we then average over translates to 
generate our hamiltonian K{s). In this case, however, 
the initial interaction term is allowed to have support 
equal to all of V. The only constraint we make is that 
k{s) must decay rapidly which means that fc(s) can be 
written as a sum; 

m — 1 

k{s) = E k^{s) (4) 

where supp(fcQ,(s)) = A^, and A^ consists of all the sites 
within a distance a of site 0, i.e., A^ = {j \ d{0,i) < 
a}. As a result, ka{s) is an operator with a support (or 
"radius") consisting of a sites centred on site 0. The 
rapid decay condition is then that 

||fca(s)|| < /(a), 0<a<m. (5) 

where f{a) is some rapidly decreasing function of a. 

We say that a hamiltonian K{s) constructed from the 
interaction fc(s) has rapid decay. We write the final 
hamiltonian resulting from this construction as 

'm— 1 

^(^) = EE%"(^)' (6) 

where fcj,„(s) = T^^" (Ti^ (fc„(s))). 

Finally, we set out the problem we aim to solve. We 
suppose H[s) is a strictly local parameter-dependent 
hamiltonian for a 2D lattice of the form with inter- 
action h{s) having 0(1) norm and 0(1) support. We as- 
sume that the ground state |f^(s)) is unique and, further, 
that the spectral gap AE{s) between the ground state 
and first excited state satisfies AE{s) > A, Vs G [0, 1], 
where A is a constant independent of n. Finally, we 
suppose that expectation values of arbitrary local oper- 
ators A G Al, with 0(1) support, in the initial ground 
state |fi(0)) can be computed efficiently, i.e., ujo{A) — 
(r2(0)|A|r2(0)) can be computed efficiently for all A G Al- 
This would be the case when, for example, i?(0) is any 
regular classical hamiltonian, that is, [/ij(s), /ik(s)] = 0, 
Vj,k G L. Alternatively, this occurs when H{s) has a 
ground state which is exactly a 2D finitely correlated 
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state. (When H{s) has spatially varying interactions we 
must require that the ground state of H{0) is a known 
product state. We need to do this in order to avoid 
the constructions of [sS] , which show that computing the 
ground state of a disordered classical systems is at least 
NP-hard.) 

Our approximation problem is therefore the follow- 
ing. First fix some error e. Then our problem is to 
find an efficient computational method to compute, for 
any local operator A with bounded support [43|, uni- 
form approximants t-^'siA) to the exact expectation val- 
ues uJs{A) = {i}{s)\A\V,{s)) . That is, our problem is to 
efficiently compute uj's{A) so that —^s{A)\ < e for 

all s G [0, f ] and for all bounded local operators A with 
bounded support. 

The constraint that the observables whose expectation 
values are to be simulated must have bounded support 
stems from the condition that in the large-n limit such 
operators should be elements of the quasi-local algebra 
Al- We lose no generality in this assumption when ap- 
plying it to the simulation of quantum algorithms be- 
cause the answer that the algorithm computes should be 
encoded in the ground state in such a way that it can 
be read out from the expectation value of a local opera- 
tor, ft is also worth noting that any correlation function 
involving a bounded number of subsystems satisfies our 
definition of having bounded support. 

Before we end this section we introduce some notation 
for approximations. Because we have occasion to refer to 
functions for which only bounds on growth, derivatives, 
etc. are known it is convenient to adopt the following 
notation. If we have two quantities A and B then we 
use the notation A < B to denote the estimate A < CB 
for some constant C which only depends on unimportant 
quantities. In almost all the cases we consider the only 
important quantity is n, the total number of spins. Thus, 
unless we indicate otherwise, A < B means that A < CB 
for some C independent of n. Because we'll be interested 
in the consequences of allowing the minimum gap A to 
depend on n we'll explicitly retain any dependence on A 
in our calculations. 



III. EFFECTIVE LOCAL DYNAMICS FOR 
EXACT ADIABATIC EVOLUTION 

In this section we study exact adiabatic evolution for 
quantum spin systems. We show that if there is a gap 
throughout the evolution then the exact adiabatic evo- 
lution is equivalent to unitary dynamics generated by an 
approximately local hamiltonian. 

We consider adiabatic quantum evolution generated by 
H(s) as s is varied adiabatically from s = to s = 1. 



Thus we would like to understand the ground state |f^(s)) 
of H{s). We do this by setting up a differential equation 
for \n{s)}: 



d 



^\m)^p'{sms)): (7) 

where P'{s) = ■£tlr2(s))(f^(s)|) and we've set phases [Hi 
so that (f2'(s)|J7(s)) = 0. Because P'{s) is not antiher- 
mitian the dynamics generated by this equation are not 
unitary. 

There are at least two ways to set up differential equa- 
tions for |fi(s)) which do generate unitary dynamics. The 
first is via exact adiabatic evolution (see |33l . |33 | for a 
rigourous discussion of rather general results about ex- 
act adiabatic evolution): 



ds 



\n{s)) = -[p{s),p'{s)\\m)- 



(8) 



Because of the gap condition on H{s), the "hamiltonian" 
[P{s), P'{s)] for this dynamics is given by first-order sta- 
tionary perturbation theory: 



[p{s),p\s)] = \n{s)){nis)\ 



dH{s) 



ds n{s)I - H{s) 
dH{s) . 



n{s)I - H{s) ds 



-\nis)){nis)\, (9) 



where fl{s) is the ground-state energy of H{s), and 
we define „/ ^ via the Moore-Penrose inverse: 

£l(s)K— ii (s) 

The other way, which we call effectively local exact adi- 
abatic evolution, is obtained by rewriting P'{s). We ex- 
ploit the fact that H{s) has a spectral gap to find 



Pis) 



X-,{t)e-''^^'^e''"^'Ut, 



(10) 



where X-yi'^) is an even real function whose fourier trans- 
form x-y is , has compact support in [—7,7], and is 
normalised so that £7(0) = 1. (See Appendix \K\ for a 
description of C°° cutoff functions and their properties.) 
We must set 7 < A to ensure that only the ground state 
appears on the RHS of The formula ^ for P{s) 

may be verified by writing e'*^ in its eigenbasis and ex- 
ploiting the L2 unitarity of the fourier transform. 
We next use the Duhamel formula 



ds 

to rewrite ([7]): 



ds 



du, 



J 



dn{s) 

ds 



txyit)dt\n{s)) + i / Xyit)e 



His 



dHjs) 
ds 



du]e'*"^'Ut\n{s)), (11) 
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where Tu^'\m) = e™^(«)Me-™-f^('*) . Using the fact 
that x-y {t) is an even function of t and canceUing phases 
we obtain 

(12) 

By integrating this expression for ^|r2(s)) in the energy 
eigenbasis of H[s) and using the assumed gap structure 
one can find that this expression is equivalent to the usual 
expression obtained from first-order perturbation theory: 



(13) 



Thanks to our assumed form of H{s) = Hq + sH', 
with H' = Y^ieL fT-'j = EjeL '^v" iV^ C*')), we notice that 



dHjs) _ 



J2\eL^'\^ and we write 



The way to see that K{s) is approximately local is to 
use the standard Lieb- Robinson bound [IB, 113, [H, HI] . 
The Lieb-Robinson bound reads 



\A),B]\\ < |r|e-'"i(^^n(e«l*l _ 1), 



(15) 



for any two norm- 1 operators A G Ax and B G Ay , with 
{x} n F = which are initially separated by a distance 
d{x, Y). The constants v and k are independent of n and 
depend only on ||/i(s)||, which is an 0(1) constant. 
What we do is define 



fco(s)=^>''\/i') 



(16) 



and 



fc,(s) = T!'^-^'\h')^Tf''-'^"\h'), 0<a<m, 



where we define 



(17) 



-\n{.s))=zJ2T,ihims)), (14) 

with initial condition that |f^(0)) is the 
ground state of H{Q) and where Ts{M) = 

The equation tells us that |ri(s)) can be ob- 

tained from |f2(0)) by unitary dynamics according to 
the time-dependent hermitian hamiltonian K{s) — 
T,jeL^s{h'^ = EjGL^j(s)> where we write kj{s) = 
J-'s(/ij). We also write k{s) = J-s{h') for the interaction 
term k{s) which generates K{s). Furthermore, we claim 
that K{s) is approximately local for all s G [0, 1]. 



with 



Xf(t) [ I T^^-^'\M)dujdt, 

(18) 



(19) 



where — {j|(i(0,j) < a}. Obviously ka{s) has 
support supp(fca(s)) = Aq, + supp(/i'). Also note that 
k{s) — X]™=o^ ka{s) (recall that m is the diameter of the 
lattice). 

We now show how the Lieb-Robinson bound provides 
an estimate on the decay of ||A:a(s)||. Firstly, we rewrite 
the Lieb-Robinson bound (flSl) so that it is more useful: 



Ha, 



Ha 



{A)~r, ''-'iA)\\^ 



f dt' r"'-^ ([i/A„ - H^^_, , rfji? (A)]) 
Jo 

< /'*'di'||[i7A„-iJA_„r,>(A)]|| 

^0 

r\t\ 

<2 dt' \\[Ha^ - Ha^_, 
Jo 



(20) 



^ — va+K\t'\ 



where we used the fundamental theorem of calculus to 
get the first line, the triangle inequality and unitary in- 
variance of the norm to get the third line, we substituted 
the Lieb-Robinson bound (jlSp in fourth line, and we in- 
tegrated the bound to get the fourth line. The a term 
in the fourth line comes from the fact that the operator 



— -ffAc-i consists of a terms (the number of terms 
crossing the boundary). The Lieb-Robinson bound, in 
this form, says that the evolution of A with respect to 

is almost the same as that for H\^_-^ , i.e., the bound- 
ary effects are unimportant for short times. 
Now consider 
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< 2 / 

Jo \Jo 

<2/ IX7WI 

CQ 



) dt 



i{2\\h'l, 



}du I 



(21) 



< 



< a j e" 
'0 



;u|/-i 



> 1 



where to get the second Unc we apphed the triangle 
inequality, to get the third line we applied the Lieb- 
Robinson bound in the form (|20p with a — diam(AQ) + 
const, (we've dropped the dependance of the interactions 
(s) on the parameter s because for these inequali- 
ties the evolution is independent of the parameter s), in 
the third line we've broken the integral into two pieces 
and applied the different regimes of the Lieb-Robinson 
bound separately with c some constant to be chosen 
later, and in the final line we applied the decay estimates 
on X7 {t) (see Appendix for a derivation of these esti- 
mates). Thus, by choosing c < v/k we see that ||fca('5)|| 
is decaying faster than the inverse of any polynomial in 
a for a > I/7, i.e., for a > c/A, where c is some con- 
stant. In this way we see that exact adiabatic evolution 
can be thought of as unitary dynamics according to the 
paramater-dependent hamiltonian K{s) which is approx- 
imately local with respect to the metric d on the lattice. 
For an illustration of the interactions of K{s) see Fig- 
ure [TJ 

IV. EFFICIENT SIMULATION OF ADIABATIC 
EVOLUTION 

In this section we apply a Lieb-Robinson bound to 
show that dynamics according to effectively local exact 
adiabatic evolution keep local operators approximately 
local and hence show that expectation values of local 
operators in adiabatically evolved ground states can be 
computed efficiently. 

Recall that we can write the ground state |ri(s)) by 
integrating ([14]) as 

\n{s))=U{s;OmO)), (22) 

where 

Uis]0) = Te'^oKis')ds' ^ (23) 



and T denotes the time-ordering operation. Our objec- 
tive is to uniformly approximate 

iUs{A) = (17(0)|Wt(s; 0)AU{s; 0)\n{0)), (24) 

for all s G [0, 1]. The way we do this is to show that the 
operator A{s) = U^s; 0)AL{{s; 0) remains approximately 
local for all s G [0, 1] and use the assumed fact that llJq{B) 
can be computed efficiently for all local operators B. For 
simplicity we assume that the operator A is located at 
the origin and has support | supp(A)| = 1. It is easy to 
extend the results of this section to apply to operators 
with bounded support on disconnected regions, such as 
correlators. 

We now study the locality of A{s). What we do is 
first show that A(s) can be uniformly approximated in 
operator norm by the series of approximants 

A„(s) = v1Js;0)AVaJs;0), (25) 

where Va^(s;0) satisfies the differential equation 

(26) 

with Va„(0;0) = I and where Ka^s) = J2jeA^^sih'^ 
and Ac = {j | d{0,]) < a}. In words: the approximation 
Aa{s) is that operator obtained by evolving A with re- 
spect only to those interaction terms in K{s) whose cen- 
tres are within a distance a of A. Naturally this means 
that Ajn-i{s) = A{s). We use a Lieb-Robinson bound to 
show that \\A{s) — Aa{s)\\ is rapidly decaying. 

To show this we prove ||t^o^'*-' {A)~T^g°' {A)\\ is small 
for \s\ < 1 and large constant a where T^}f\M) = 
h{^{s; s')MU{s\ s'). To make this expression easier to deal 
with, and to more explicitly relate it to group-velocity 
bounds, we rewrite it: 
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FIG. 1; Illustration of the rapidly decaying interactions for the effectively local hamiltonian for exact adiabatic evolution. 



\\r::^'\A)-r::t^'\A)\\ 



(27) 



< 



ds'\\[K^As').r^}>\A)]\\, 



where Ka^o{s) = EjeL\A„ ki^)- 

We now apply a general Lieb-Robinson bound recently 
proved in llGll . In order to apply the Lieb-Robinson 
bound of [l6| we need to establish that our hamiltonian 
K{s) satisfies the conditions of Assumption 2.2 of p^ . 
which in our case reads 

m—1 

\\k4{l + 2a{a + l)f{2 + 2ar<si, (28) 

q:=0 

where 77 is a positive constant and si is some constant. 
We need to ensure that the sum on the left evaluates 
to a constant instead of diverging. The only flexibility 
we have is to choose a decay estimate for which 
is strong enough to overwhelm the polynomial in a it is 
multiplied by. The highest power of a appearing in this 
sum is a^^^ . Thus we use the decay estimate (PT|) and 
choose I > 7 + r], and so we find that this constant si 



equates to 

si = a(0/y, (29) 

where 

'^^ (l + 2a(a+l))^(2 + 2a)-' 
<^{l) =ci2^ ^TTT , (30) 

a— 1 

and c/ is the constant arising from the estimate (|2ip . and 
I is any chosen power I > 7 + rj. (The proof of the gen- 
eral Lieb-Robinson bound described in is easily ex- 
tended to cover parameter-dependent hamiltonians such 
as K{s).) This reads 

p{l)\Y\ (e^l^l-l) 

(31) 
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for any two norm-1 operators A G Ax and B G Ay, with 
{x} nY — $ which are initially separated by a distance 
d{x, Y) and p{l) is a constant which depends only on 
The constant v is independent of n and depends only on 
||/i(s)l|. We isolate the dependence of this bound on the 
minimum gap 7 by defining 



(32) 



Note that we are going to systematically redefine this 
function in our subsequent derivations to absorb extra 
constants and occurrences of 7. With this bound we first 
obtain an upper bound on \\[TfJf\A),ky„{s')]\\ (recah 
that the operators kj,a{s) are defined via Eq. fH)): 

\\[A{s),ki4s)]\\< I (1+5-")' 

(33) 

where S = d(0,j) and (1 + 2a{a + 1)) = |Aa(j)| = 
{x I (i(j, x) < a}\. We use the estimate (|21l) and redefine 
5(7, 1) to find the upper bound 



• g(7,0(l+2a(a+l)) 



< 6 



- \2||A||l|fc„(s)||, a>6. ^ ^ 

We next find the minimum of the denominator a''^'^{l + 
S — aY on the interval 1 < a < 5, which is 5', and redefine 
g{'y,l) to arrive at the final upper bound 

||[A(s),fcj,„(,s)]||<|3£|'l (35) 



Thus, by choosing the centre j far enough away from 
the centre of A{s) we find the behaviour 



\\[A{s),kUs)]\\ 



< 



.9(7,0 
d(OJ)' 



yi > 1, 



(36) 



i.e., the quantity ||[A(s), A:j_ct(s)]|| decays faster than any 
polynomial in d{0,j). 

We next use our upper bound psp to obtain an upper 
bound on ||[A(s), A:j(s)]||: 

m — 1 

\\[A{s),ki{s)]\\<J2\\[A{s),kUm 



< iM , V (37) 



7'T"-^a' 



< 



g(7,0 
§1-1 ' 



where we've redefined 5(7,/) in the last line. 

Now we use the decay estimate ([57]) in (|27p to provide 



an upper bound for \\A{s) — Aq(s)||: 

\\Ais)-A^{s)\\< [\s\\[A{s),k^is)]\\ 

""^ (l + 2^(^ + l))g(7,0 (38) 



(5 — Q 

5(7,0 



< 



where we've redefined 17(7, 1). 

So, as long as a is chosen to be so large that it over- 
whelms the 0(1) constant g{'y,l) we find that \\A{s) — 
Aa{s)\\ can be made to decay faster than any polyno- 
mial in a, and hence, can be made as small as de- 
sired. Thus there exists some constant a such that 
||A(s) — Aa(s)|| < e. Note that, because k{s) has support 
throughout L, Aa{s) has support throughout L. 

In order to provide a simulation method to compute 
approximations to ground-state expectation values (A) 
we need to show that Aa{s) can be approximated by an 
operator with support only on a constant number of sites 
around supp(A) — 0. The way we do this is to show that 
Aa (s) is operator- norm close to 

A^As)^vIJs)AVa^Js), (39) 
where Va^ „ satisfies the differential equation 



' E •^'l''"'('^j)VA„..(^) = ^KA^As)VA^Js), (40) 



and 



(41) 

with A/3(j) = {x|ji(j,x) </3}. 

To show that Aq,_^(s) is close to Aa{s) we first exploit 
the general inequality 

I|VaJ.)-Va„,,(s)|| < t\\KAAs') - KA^As')\\ds' 
Jo 

(42) 

which is proved, for example, by exploiting the Lie- 
Trotter expansion, and then upper-bound the right-hand 
side using the triangle inequality by 



\\KA^{s')-KA^,{s')\\ds' < 



E t\\his')-kAs')\\ds\ (43) 
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where fcj,^(s) = T-^^ "'"(/i^). We can upper-bound the 
integral on the right-hand side by using an argument 
identical to the one used to show pTjl . We thus obtain 

jeA. jeA„ ^ ^ (44) 



where I is any power, and we've used the fact that the 
number of sites in is given by 1 -I- 2a{a + 1). By 
choosing > a we find that Va^ (s) can be made as close 
as desired to Va„ ^ (s)- 

To obtain closeness of our final approximation Aq,^^(s) 
to A{s) we use the triangle inequality 

\\A{s) - < \\A{s) - A^{s)\\ + \\A^{s) ~ A^^p{s)\\ 

(45) 

where we've used the upper bound (|38p with an adjusted 
value of I and we've also used ((44|) with an appropriate 
choice of power I' . We therefore find that it is sufhcient, 
for a given constant e to choose large (but 0(1)) a and 
/3 so that 

\\A{s) - A^{s)\\ < e. (46) 

The actual values of a and /? required to reduce the 
error to below e scales better than linearly with 
w — max((7(7,^), I/7), where 7 is a constant multiplied 
by the minimum energy AE encountered along the adia- 
batic path. Thus the support of the final approximation 
Aa,/3is) is given, in the worst case, by supp(^Q,^_a(s)) < w. 
Note that w depends, via 3(7, Oj exponentially on I/7, 
i.e., the inverse energy gap. 

Because the final approximation Aa^p^s) can be com- 
puted via integrating (gD]), and by noticing that this 
integration can be performed by restricting our atten- 
tion to the finite-dimensional subalgebra Aw, where 
W = supp(y4Q_^(s)), we see that Aq_^(s) can be com- 
puted using resources which scale as 2^^™, with c some 
constant. 



V. DISCUSSION 

In this paper we have shown how to efhciently cal- 
culate the ground-state expectation values of local op- 
erators with constant support for gapped adiabatically 
evolving spin systems. In order to provide our simulation 
method we reduced the problem to showing that under 
exact adiabatic evolution the expectation value of a local 
operator can be computed from the expectation value of 
an approximately local operator in the unevolved ground 
state. Given this observation we then argued that if it 



is easy to compute expectation values of local operators 
in the original ground state then one could approximate 
the desired expectation values arbitrarily well by using 
time and space resources that scale with the inverse gap. 

Our approach has several shortcomings. The first is 
that the scaling of the simulation resources with the error 
e scales faster than 2^/*^. This means that if the expec- 
tation value of an operator which is a sum of many local 
operators is desired then our simulation method may re- 
quire superpolynomial resources. For example, if the ex- 
pectation value of the total magnetisation M = Sjgl '^j 
(as opposed to the more traditional average magneti- 
sation m = M/n) is required to some accuracy e then 
our simulation method will require superpolynomial re- 
sources. This is not entirely unexpected, after all, in the 
thermodynamic limit such operators are unbounded and 
cannot be approximated at all. Another manifestation of 
this shortcoming is that if the expectation values of the 
local operators are required to an accuracy which scales 
as e < 1/n then our method may require superpolyno- 
mial resources. These problems do not manifest them- 
selves for the applications we have in mind. Namely, 
when applied to the calculation of average properties of 
two states in the same quantum phase we only require ac- 
curacy to some small constant e which doesn't scale with 
the system size, and when applied to simulating adiabatic 
quantum algorithms we only need e to scale as a constant 
in order to read out the answer of the algorithm. 

The second shortcoming of our method is that, by the 
current method, we are unable to directly approximate 
the scaling of the geometric entropy [s^l 'S'a with A. The 
reason for this is that our current method approximates 
Pa(s) by calculating approximations to all the expecta- 
tion values of a basis of operators for ^a- Because we 
are computing approximations to expectation values we 
end up computing only an approximation pi\{s) to pa{s). 
The best continuity result available for the von Neumann 
entropy is Fannes inequality (see, for example, [2^) for 
a derivation) which implies that the error in the approx- 
imation S'a calculated from pa(s) grows larger as A in- 
creases. We'll describe an approach to this problem using 
exact adiabatic evolution in a future paper. 

The principle characteristic of our approach is that ap- 
proximations are made in the Heisenberg picture. What 
we mean here is that instead of approximating the 
evolved quantum state of the spin system in operator 
norm we instead compute approximations to the evolved 
local operators. We should expect this strategy to be 
successful because the locality of the interactions in the 
hamiltonian doesn't manifest itself in the Schrodinger 
picture but, thanks to the Lieb- Robinson bound, it is pre- 
cisely clear what locality implies for local operators in the 
Heisenberg picture. Because in the thermodynamic limit 
we are only able to physically access local operators (such 
as average magnetisation and correlators) this approach 
doesn't lead to any loss of generality over computations 
carried out in the Schrodinger picture. 

It is possible that our analysis acutally applies to all 
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Phase boundary = quantum 
phase transition; poly(«) 
storage resources required. 
Time complexity? 



Mystery zone (conjectured); 
requires at least subexponential 
resources? 




Classical spin model point; 
poly(«) resources required 



Yarotsky boundary 



FIG. 2: Conjectured diagram of the space of local spin models and the associated computational resources required to compute 
approximations to local ground-state properties. 



gapped spin models. This is because it is possible that 
any gapped spin model is adiabatically connected [sij to 
a classical spin model with trivial ground state. Classi- 
cal renormalisation-group style argumentation certainly 
seems to back this statement up: after all, we know that 
the RG fixed points are either trivial (classical) or quan- 
tum critical points. However, there is as yet no rigourous 
general proof of this statement for quantum spin systems. 

We would like to suggest that the following description 
of the space of local (translation-invariant) spin models 
is correct. Firstly, in this space there are many distin- 
guished points, classical spin systems, where the ground 
state can be calculated trivially. Around each of these 
points is a small region in hamiltonian space of hamilto- 
nians which are provably adiabatically connected to the 
classical spin model points 0, [H, [1^ . In these regions 
we have shown that the local ground-state properties can 
be determined efficiently. Outside these small regions 
there are other regions which may or may not be adi- 
abatically connected to the classical spin model points 
where the hamiltonians are gapped. In these regions it 
is known that the local ground-state properties can be 
calculated using subexponential resources [l3|. On the 
boundaries between the quantum phases there are quan- 
tum critical walls. For these points, in ID, it is known 
that an approximation to the ground state as a finitely 
correlated state can be stored using polynomial space 
p^ . It is not known if these approximations can be ob- 



tained efficiently. This picture is summarised in Figure [H 
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APPENDIX A: PROPERTIES OF SMOOTH 
CUTOFF FUNCTIONS 

In this Appendix we briefiy review the properties of 
compactly supported C°° cutolf functions. 

Of fundamental utility in our derivations is a class of 
functions known as compactly supported C°° bump func- 
tions. These functions are defined so that their fourier 
transform x-yi'-^) is compactly supported on the interval 
[—7, 7] , and equal to 1 on the middle third of the interval. 
Such functions satisfy the following derivative bounds 



(Al) 



for all j with the implicit constant depending on j. This 
is just about the best estimate possible given Taylor's 
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theorem with remainder and the constraints that x-yi^) 
is equal to 1 at w = and X-yi^) is compactly supported. 

The function Xy {t) has support throughout K but it is 
decaying rapidly. To see this consider 

X.W = -^£^e-A^,H<i. (A2) 

which comes from integrating by parts. Continuing is 
this fashion allows us to arrive at 



(|Aip . and using the compact support of ^^(i^) we find 



\Xt{t)\ 



< 



7 



r7 1 
< / ,duj 



< 



\lt\' 

1 



(A4) 



X7W 



1 

2^ 



Xj{uj)(hj (A3) 



Since x~f{uj) has all its derivatives bounded, according to 



for all i G N. Thus we find that x-y (t) decays to faster 
than the inverse of any polynomial in t with characteristic 
"width" 1/7. The existence and construction of such 
functions is discussed, for example, in [s^ [37j. 
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While the analysis of was carried out for ID systems, 
it is clear how the calculations generalise to 2D systems 
and show a subexponential resource scaling to obtain a 
representation of the ground state as a 2D finitely cor- 
related state fXL, ^] or "projected-pair entangled state^' 
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lem to determine the classes of 2D PEPS for which it 
is computationally easy to extract local properties, such 
as correlators. Indeed, because of certain classical con- 
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At least those which are equivalent to conformal models. 
Note that it is easy to construct models which are "crit- 
ical" in the sense of vanishing spectral gap, but which 
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be decided in polynomial time on a quantum computer. 

[45] This idea has been recently explored in ^39,] where it was 
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quantum algorithms. Note that the codes investigated in 
this paper do not increase the overall gap of the original 
hamiltonian after the encoding. 

[46] The order notation 0( ) refers to the growth or decay of 
quantities with respect to the fundamental parameter of 
this paper: n, the number of spins. 

[47] That is, we assume |supp(^)| is any (arbitrarily large) 
0(1) constant. 

[48] Note that this choice precludes an extension of our anal- 
ysis to study Berry phases in spin systems. We'll relax 
this constraint in a future paper. 

[49] Because of our < notation we'll reuse the symbol c by 
systematically redefining it. 

[50] Recall that the geometric entropy Sa is equal to the 
von Neumann entropy of the restriction of the ground 
state p(s) = |n(s))(J7(s)| to a contiguous block of A 
spins: Sa = S'(pa(s)), where S{p) — — tr{p\og2{p)) and 
Pa{s) — tr^(/9(s)), with tr^ denoting partial trace over 
all spins except those in A. 

[51] What we mean here is that there exists some adiabatic 
path of local hamiltonians from a local classical spin 
model to the desired gapped spin model where there is a 
constant gap along the entire path. 



